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We present results of numerical simulations of coupled Ginzburg-Landau equations that describe 
parametrically excited waves. In one dimension we focus on a new regime in which the Eckhaus 
sideband instability does not lead to an overall change in the wavelength via the occurrence of a 
single phase slip but instead leads to "double phase slips". They are characterized by the phase 
slips occurring in sequential pairs, with the second phase slip quickly following and negating the 
first. The resulting dynamics range from transient excursions from a fixed point resembling those 
seen in excitable media, to periodic solutions of varying complexity and chaotic solutions. In larger 
■ systems we find in addition localized spatio-temporal chaos, where the solution consists of a chaotic 

region with quiescent regions on each side. We explain the localization using an effective phase 
diffusion equation which can be viewed as arising from a homogenization of the chaotic state. In 
• two dimensions the double phase slips are replaced by fluctuating bound defect pairs. 

Q 

\Q ■ I. INTRODUCTION 

Pattern formation is a striking and interesting phenomenon observed in many physical systems driven far from 
. thermodynamic equilibrium. In addition to parametrically excited waves (which are described below), classic examples 
include the rolls formed in Rayleigh-Benard convection when a fluid confined between two horizontal plates is heated 
Cn| , from below, the annular Taylor vortices that form in a fluid between rotating concentric cylinders, and rotating spirals 
(and other patterns) consisting of regions of different chemical concentrations formed in chemical reactions such as 
the Belousov-Zhabotinsky reaction. A review of pattern forming systems can be found in Typically, the basic 
equations that describe these pattern- forming experiments are quite complicated (such as the Navier Stokes equations). 
However, for pattern-forming systems near onset, the complicated system of equations can often be simplified using 
asymptotic expansions for small amplitudes. The resulting amplitude equations (e.g. ^-^]) describe the variation of 
the pattern on large space and time scales. 

In many systems transitions from simple to chaotic behavior are found when a control parameter (such as the 
temperature difference in Rayleigh-Benard convection, the angular velocity of the inner cylinder in Taylor-Couette 
flow, or the supplied concentration of a particular chemical in a reaction diffusion system) is changed. Studies of 
chaos began with low-dimensional attractors as found, for instance, in the Lorentz system |(| and have progressed to 
chaotic behavior in both space and time in spatially extended systems with a large number of degrees of freedom. Of 
particular interest in both experimental and mathematical systems is identifying features that allow the different types 
of observed chaos to be classified in some way. A well studied example is the complex Ginzburg-Landau equation, 
for which a number of different chaotic regimes have been identified 0-0. Relevant for the present paper is the 
observation that a distinction can be drawn between a phase-chaotic regime and amplitude-chaotic ones. In one 
dimension amplitude chaos is characterized by the occurrence of phase slips during which the amplitude of the wave 
goes to zero and the total phase of the system changes by 2ir, that is, a wavelength is inserted or eliminated. In 
phase chaos there are essentially no phase slips |9|,|llJ and the system can be described by an equation for the phase 
of the wave alone |l^ ] . For the complex Ginzburg-Landau equation this phase equation is the Kuramoto-Sivashinsky 
equation. Phase equations break down during phase slips so they cannot be used to describe amplitude chaos. 

In this paper we consider parametrically excited waves, which are also a classic pattern forming system [ p"3[ . 
Parametrically excited waves arise quite generally in systems exhibiting weakly damped modes that are oscillatory in 
both space and time, when these modes are forced at twice their natural frequency. Faraday's experiment []l3f in which 
surface waves are formed on a liquid in a vertically oscillated container is probably the most well known example. 
Stripe, square, and hexagon patterns (e.g. |[4|), as well as quasipatterns pHf , spirals and chaotic states (e.g. 
P?] , ^8|jl4| ) have all been observed in Faraday experiments with various fluids, forcing amplitudes, and frequencies. 
Parametrically excited waves have also been observed in magnetic materials (spin waves) [|l9| , po| , in electroconvection 
in nematic liquid crystals (2^j2^| , and in Taylor-Dean flow [^3) . 

In this paper we discuss dynamics arising in numerical simulations of a system of equations that model parametrically 
excited waves and focus on solutions which are chaotic in space and time. In our simulations we find chaotic solutions 
which are characterized by a new phenomenon: double phase slips in one dimension and the analogous bound defect 
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pairs in two dimensions. These phenomena are described in sections III and VI. In addition to chaotic solutions where 
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these phenomena occur irregularly in space and time, which are describ ed in secti ons [IV C| and V A, we also describe 
transient and periodic solutions involving double phase slips in sections [V A and IV 

Of central interest in our one-dimensional simulations are localized chaotic states (section VB). That is, there 
are solutions that include a spatial region, bounded by quiescent regions on each side, in which the dynamics are 
irregular in space and time. Experimentally, similar phenomena have been observed in Taylor vortex flow pi] , 
Rayleigh-Benard convection |2^], and parametrically excited surface waves [ pj|]l4| l. The localization mechanism in 
these experimental systems is only poorly understood. For the system of equations discussed in this paper the 
localization can be understood using an effective phase diffusion equation. We confirm the validity of this approach 
by explicitly determining the effective phase diffusion coefficient in section |V C| . A short account of these results has 
been presented previously [gfjj. 

While most of the results in this pape r are for a one dimensional system, section [Vl| describes results for an analogous 
system in two dimensions. Section VII concludes our paper with a summary and discussion of the results. 



II. THE MODEL 



The equations we consider in this paper generically describe a system with a uniform state that has a supercritical 
Hopf bifurcation to a state that is periodic in both space and time when subjected to a uniform, constant forcing. 
Subjecting the system to an additional spatially uniform, time-periodic forcing at approximately twice the Hopf 
frequency creates a codimension-three bifurcation point involving the following three parameters: 1) the amplitude of 
the time-independent forcing, 2) the amplitude of the periodic forcing, and 3) the frequency of the periodic forcing. 
In the vicinity of this codimension-three bifurcation point the following amplitude equations |p7[-p9| are derived by 
considering the reflection and translation symmetries in space and the weakly broken translation symmetry in time: 

d T A + sd x A = dd%A + aA + bB + cA(\A\ 2 + \B\ 2 ) + gA\B\ 2 , (1) 
d T B - sd x B = d*d\B + a*B + bA + c*B(\A\ 2 + \B\ 2 ) + g*B\A\ 2 . (2) 

The dependent variables A and B are complex, and represent amplitudes of left and right traveling waves that are 
summed together to yield the solution of the underlying system: 

u(x, t) = eA(X, T)e t{q " x "^ t] + eB(X, T)e i{ - q ' x+ ^ + c.c. + o(e), (3) 

where u(x,t) would give, for example, the position of the fluid surface in the Faraday experiment. The amplitudes 
A and B vary on the slow time and space scales, T — e 2 t and X — ex. The coefficients in (0J||) are complex except 
for s and 6, which are real. The real part of the coefficient of the linear term a gives the linear damping of the 
traveling waves in the absence of the periodic forcing and is proportional to the distance from the Hopf bifurcation. 
The coefficient of the linear coupling term b gives the amplitude of the periodic forcing. (Note that this term breaks 
the time translation symmetry.) The imaginary part of the coefficient of the linear term a gives the difference between 
the frequency of the unforced waves and half the forcing frequency uj e . 

While in general a codimension-three bifurcation might be difficult to realize in an experiment, in the above situation 
(involving a system that undergoes a Hopf bifurcation) each of the three bifurcation parameters are natural control 
parameters in the system. In the Faraday experiment, however, there is no Hopf bifurcation and the above equations 
are only valid in the limit of weak damping (small viscosity). In this case the coefficients of the nonlinear terms c and 
g would be purely imaginary, since their real parts also represent (nonlinear) damping and would be higher order. In 
addition, for the Faraday experiment the group velocity parameter s would generally be complex, and s in (|^) would 
be replaced by its complex conjugate s*. 

There is one caveat concerning when they are used to describe either the Faraday experiment or parametrically 
excited waves induced near a Hopf bifurcation point. With the scaling X — ex, the group velocity term (involving 
dx) would generically appear at lower order, so (0J|) require that the group velocity is small. If it is not small a strict 
asymptotic analysis would require the introduction of a second slow time scale and would lead to equations in which 
the cross-coupling between A and B arises through a spatially averaged term fl30|-|34t . 

Note that (|l]|2|) are in one (spatial) dimension. With periodic boundary conditions (as will be imposed in our 
numerical simulations) these equations would describe (for example) the Faraday experiment in an annular container. 
Equations in two (spatial) dimensions (see (|27|j28| ) below) are discussed in section VI . 
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Solutions of (0,^|) of the form A — Aoe*"® , B = Boe l< ^ x include those where Aq and Bq are constants of unequal 
magnitude that correspond to traveling waves in the underlying system (see (j^)) , and those where Aq and Bq are 
constants with equal magnitude that correspond to standing waves in the underlying system |2?],|28|]. It is these 
standing-wave solutions, which are phase locked to the periodic forcing, that are the focus of this paper. In particular, 
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we discuss the ensuing dynamics when these standing waves are perturbed. There are also solutions where Aq and Bq 
have equal magnitude but vary periodically in time, which correspond to standing waves that are not phase locked 
to the external forcing p7[ . We do not treat these standing waves in this paper. 

The standing waves which are phase locked to the external forcing are stationary solutions of (0J^) and are given 

by fH 

A = r^(Q x +M ( 4 ) 

B = ife*« x +*») (5) 
2 _ -(a r n r + ajUi) ± ^{a r n r + a l n i ) 2 - (n? + nf)(a% + a 2 - b 2 ) 

6: " {n 2 r + n 2 ) { > 
i i , h + n.R 2 

9b- 4>a = arctan - — ■ — (7) 

a r + n r H z 

and: a = a r + ia, = a — isQ — dQ 2 (8) 

n = n r + Wij = 2c + g. (9) 

They arise from the trivial solution in a steady bifurcation at the onset b 2 — a 2 + a 2 . Equations (^,||) actually give a 
family of solutions for each wavenumber Q since the phase 4>a (or 4>b) can be chosen arbitrarily. Changing the value 
of 4>a (with 4>b — 4>a fixed) corresponds to spatial translations of the solution. 

Close to onset (fflfl) can be reduced to a single Ginzburg-Landau equation with real coefficients |3j| . This implies 
that the waves are stable near onset in a band of wavenumbers which is limited by the Eckhaus (i.e. a long-wavelength) 
instability. Near onset the neutral and Eckhaus curves are parabolic [[36| . For larger values of the forcing amplitudes 
b, the Eckhaus curve need no longer be parabolic and can assume a variety of shapes. For a — —0.05, c = — 1 + 4i, 
d = 1 + 0.5i, s = 0.2, g = — 1 — 12i, the Eckhaus curve closes on itself as shown in Figure | @. Thus, for b > 0.85 all 
standing waves become unstable, somewhat reminiscent of the Benjamin- Feir instability of traveling waves |37|. In 
this paper we focus on the dynamics resulting from the Eckhaus instability. All computations presented in this paper 
were made using the same parameter values as used in Figure [l]. 

In reference [ p5| the Eckhaus curve in Figure [l] was obtained using phase dynamics. A more complete linear stability 
analysis considering sideband perturbations with arbitrary wavenumbers confirms that there is no short-wavelength 
instability for these parameter values. 



III. PHASE SLIPS AND DOUBLE PHASE SLIPS 

The main questions to be addressed in this paper are: What are the solutions to (QJ^) if started with initial conditions 
that correspond to unstable standing waves? What features characterize the solutions? How do the solutions differ as 
one changes the forcing amplitude b and the wavenumber of the initial condition Q? We investigate these questions 
by calculating numerical approximations to the solutions of ([l],^]) using a pseudospectral method in space (periodic 
boundary conditions) and a fourth-order Runge-Kutta/integrating-factor scheme in time. 

The dynamics described in this paper all incorporate a feature we call a "double phase slip" . In this section we 
describe this feature and contrast it to an ordinary single phase slip, which is a well known feature of many systems 
(including the solution of a single Ginzburg-Landau equation with real coefficients) when started with an initial 
condition that is Eckhaus unstable. In order to excite the Eckhaus instability we typically choose an initial condition 
that is a wave, perturbed so that the local wavenumber varies slightly throughout the domain such that it has a 
maximum and a minimum at distinct locations. 

Since near onset (|l],||) can be reduced to a single real Ginzburg-Landau equation p5|, the solution of (|l|,^) for 
small values of b should resemble the solution of a single real Ginzburg-Landau equation. This is indeed seen in our 
simulations. For small values of b (b = 0.1, say) and an initial condition slightly to the right of the Eckhaus stable 
region in Figure ^, the solution undergoes a phase slip near the location of the initial maximum in the wavenumber 
that reduces the wave number and moves the solution into the stable region. This behavior is shown in Figure |^a. 
For larger values of b however (b = 0.6, say), the behavior of the solution is markedly different, as shown in Figure 
^b; given an initial condition slightly to the right of the Eckhaus stable region, the solution undergoes a phase slip 
that reduces the average wave number but a short time later undergoes a second phase slip at essentially the same 
location in space that restores the wave number to its original value. This double phase slip causes the solution to 
remain in the Eckhaus unstable region, thus allowing persistent dynamics. 

In our simulations a phase slip in A is always accompanied by an essentially coincident phase slip in B. Thus when 
we refer to a double phase slip in the solution of (0^) we mean that there are double phase slips in both A and B at 
approximately the same point in space and time. 
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Double phase slips occur in a variety of solutions to (|lj,0) including periodic solutions, chaotic solutions, and non- 
trivial transients. These solutions are described in the following sections. We first present results for small system 
sizes. Then we discuss large systems and show that they allow a simplified long-wave description. 



IV. DYNAMICS OF A SMALL SYSTEM 



For the simulations presented in this section the domain size L has been chosen such that the initial conditions, 
which have an average wave number in the vicinity of the Eckhaus curve, include only five wavelengths. Three 
different types of solutions are presented: 1) Transient solutions where there is only a single double phase slip, 2) 
Solutions where double phase slips occur periodically in time, and 3) Complex solutions where double phase slips 
occur irregularly in space and time. Figure | shows where examples of these solutions occur in relation to the Eckhaus 
curve that is shown in Figure [|. The initial conditions for these solutions are close to the standing wave solutions 
but perturbed such that there is a maximum in the local wave number that triggers the first double phase slip and 
determines its location. Examples of each of the solution types are described below. 



A. Transient Solutions 



In the transient solutions described here, there is only one double phase slip. After the double phase slip, the 
solution relaxes to the stationary solution corresponding to the regular standing wave state in the underlying system. 
These transient solutions resemble the response of an excitable medium; beginning with an initial condition that is 
in the vicinity of a stable fixed point, the transient solution makes an excursion far away from the fixed point before 
finally coming to rest there. This is most clearly seen in the phase-plane representation of the solution shown in Figure 
^. It shows the magnitudes of the (spatial) Fourier modes corresponding to 4 and 5 wavelengths in the domain. Each 
dot in Figure |J represents an instant in time. The square indicates the initial condition which is the perturbed 5-wave 
state. The loop in Figure ^ is traversed once, with the final state of the solution being the fixed point on the vertical 
axis (marked by a diamond) corresponding to the (stable) 5-wave solution. 

Solutions with a single double phase slip occur only where the initial condition is a perturbation of a stable solution. 
The lowest point on Figure || appears to contradict this statement since it appears on the outside of the Eckhaus curve. 
However, this is due to the fact that the Eckhaus curve represents a long-wavelength instability. Since arbitrarily 
small wavcnumbers are not present in our finite domain the true stability curve for our finite system is slightly outside 
the Eckhaus curve and all of the solutions with a single double phase slip are in fact within the stable region. At 
the stability limit an unstable saddle-point solution branches subcritically off the branch of solutions given by (^j|) 
p8| , p9[ . The perturbation required to excite a double phase slip must place the initial condition beyond this saddle 
point in phase space. As the wavenumber is moved further into the Eckhaus stable region, the bifurcated saddle point 
moves further from the stable stationary solution in phase space, and a larger perturbation is required to excite the 
system. 



B. Periodic Solutions 



In the periodic solutions, double phase slips occur at the same locations in space, periodically in time. An example 
of the simplest of these periodic solutions is shown in Figure | where the location in space and time of each double 
phase slip is indicated by a circle. The time between each double phase slip is much larger than the time between the 
individual phase slips which make up a double phase slip. This can be appreciated by looking at Figure ^ which shows 
the (spatial) average of the wavenumber plotted versus time. Each of the downward "spikes" in Figure ^ represents 
a double phase slip; the average wavenumber is reduced by the first phase slip, then restored by the second. Figure 
shows the phase-plane representation of the simple periodic solution. The loop in Figure [7] is traversed periodically, 
with the upper left portion of the loop (where the mode 5 magnitude is large and the mode 4 magnitude is small) 
being passed between double phase slips, and the lower right portion of the loop (where the mode 5 magnitude is 
smaller and the mode 4 magnitude is larger) being passed during a double phase slip. 

Although the periodic solutions occur for larger wavenumbers (for a given value of 6) than the transient solutions 
where there is only one double phase slip (see Figure ||) we find that the periodic solutions can still occur within the 
Eckhaus stable region indicating that stable stationary solutions and stable periodic solutions can coexist. 

As the wavenumber is decreased the period of the simple periodic solutions diverges as shown in Figure ||. This 
appears to be due to the unstable saddle-point solution which branches subcritically off the branch of solutions given 
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by (ffijty at the Eckhaus stability limit. Inside the Eckhaus stable region, as the average wavenumber of the initial 
condition is decreased by increasing the size of the system L, the saddle point moves away from the stable fixed 
point and towards the periodic orbit in phase space. This causes the periodic orbit to pass closer to the saddle 
point, increasing its period (cf. Figure ||). Our simulations suggest that the divergence of the period occurs when 
the periodic orbit becomes a homoclinic orbit of the saddle point. Beyond this point the periodic solutions no longer 
exist and the excitable transient solutions occur. 

Increasing the wavenumber from those which produce simple periodic solutions we find a period doubling bifurcation 
which leads to the solution shown in Figure |y. Increasing the wavenumber further results in even more complicated 
periodic (or quasiperiodic) states such as that shown in Figure |l0| where a wave of double phase slips propagate 
through the system. The exact form of the more complicated periodic solutions may depend on the value of b but 
the trend of simple periodic solutions losing stability to more complicated periodic solutions was observed for all 
investigated values of b that exhibit periodic solutions. 

C. Complex Solutions 

Figure [ll] shows a complex solution where the double phase slips occur irregularly in space and time. The complex 
solutions occur for even larger wavenumbers (for a given value of b) than the periodic solutions. If the average wave 
number of the initial condition is too large however, a single phase slip occurs, which reduces the wavenumber. The 
phase slip need not take the solution all the way into the stable band. Instead the solution can remain in the domain 
of periodic or even complex dynamics. 

V. DYNAMICS OF A LARGE SYSTEM 

In this section we present solutions for larger domains than the solutions discussed in the previous section. For 
the particular solutions discussed, approximately 50 wavelengths with a wavenumber near the Eckhaus curve fit in 
the domain whereas in the previous section only 5 did. We focus on the complex solutions where the double phase 
slips occur irregularly in space and time. Based on the irregular appearance of these solutions we believe them to be 
chaotic. Autocorrelation functions calculated using the location of the double phase slips show that these locations 
are not correlated for large distances and times, confirming the chaotic nature of the solutions in both space and 
time. We discuss two different solution types: 1) Extended spatio-temporal chaos where the double phase slips occur 
irregularly across the entire domain, and 2) Localized spatio-temporal chaos where the irregular occurrence of double 
phase slips is restricted to only part of the spatial domain. 

A. Extended Spatio-Temporal Chaos 

Figure |l2| a shows a space-time diagram similar to Figure [ll] for an extended chaotic solution in a large system, with 
each double phase slip represented by a single dot. The extended chaotic states exist over a range of wavenumbers 
with the space-time density of double phase slips increasing with wavenumber as seen in Figure [T^J. Beyond a certain 
wavenumber the chaotic state loses stability, however. A solution with such a large initial wavenumber is represented 
in Figure |l2] c where at t w 13000 and x rj 265 a single phase slip occurs, which reduces the local wavenumber. 
Subsequently there are no phase slips in this area until the higher wavenumber in the regions on cither side diffuses 
in. Eventually the wavenumber and the density of double phase slips in Figure [T^ c become homogeneous again and 
their final values are the same as those in Figure |l2|b. The transition appears to be related to the fact that with the 
large wavenumber (that is, before the single phase slip in Figure |l2|c) the time between the occurrence of double phase 
slips at a given location in space is similar to the duration of each double phase slip, i.e. the time between the two 
phase slips within a double phase slip. This is due to both the high density of double phase slips and their tendency 
to reoccur at the same location in space. 

B. Localized Spatio-Temporal Chaos 

Perhaps the most interesting solutions arising in our simulations are those featuring localized spatio-temporal chaos. 
An example is shown in Figure |l3|. As in Figure [l^, each double phase slip is represented by a single dot. We start 
with an initial condition with a maximum in the local wavenumber that triggers an initial double phase slip. This 
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initial double phase slip triggers additional double phase slips in the vicinity of the first. Originally the size of the 
region in which the double phase slips occur grows with time. However, for later time, the size of the region ceases to 
grow and the chaos is seen to be localized. Outside of the region in which the double phase slips occur, the solution 
is the time-indcpcndcnt solution which corresponds to standing waves in the underlying system. 

Figure [l4| shows the width of the chaotic region as a function of time for different initial widths bracketing the 
final width of the state shown in Figure For both initial conditions the width converges to the same final value 
indicating that in the final state the spatio-temporal chaos is indeed localized. 

C. Effective Phase Equation 

The localization mechanism for the steady and chaotic regions in Figure |l^ can be understood using an effective 
phase diffusion equation. In this section we present this analysis. Our analysis uses the phase of A(X, T) but an 
analysis based on B(X,T) would yield the same results. 

From our investigation of solutions with extended chaos we expect the chaotic region to have a larger local wavenum- 
ber than the region where the solution is time independent (see Figure |^) . Figure [l5| shows the local wavenumber of 
the time average of the solution A(X, T) (averaged from t = 2, 000, 000 to t = 2, 010, 000) for the state represented in 
Figure ^| and confirms this; the local wavenumber in the chaotic region is larger than the local wave number in the 
quiescent region. The local wavenumber for B(X,T) is similar. So the steady and chaotic regions in Figure [To] can 
be viewed as domains with differing wavenumber. 

Somewhat similar stationary states consisting of domains of differing wavenumber have been observed experimen- 
tally in convection in narrow channels f4C(| , as well as in Taylor vortex flow p3 ] and in two-dimensional optical patterns 
ffl| . The convection states can be described by the phase diffusion equation [] [ p4||3^fl^|43| ] 

d T (f>(X, T) = D(q)d 2 x <t)(X, T) + h.o.t, (10) 

which is nonlinear due to the dependence of D on q = dx4>- The higher order terms include higher derivative terms 
(e.g. Gd x (f>(X,T)) which are important when D is small or negative, or when these higher derivatives are large. 
Stable solutions of jl0| ) that consist of domains of differing wavenumber can occur if the diffusion coefficient D(q) is 
negative only over a small range of wavenumbers. In this situation, initial conditions with a uniform wavenumber in 
this unstable range will evolve to a structure consisting of domains in which the local wavenumber lies in either of 
the two adjacent ranges of stable wavenumbers. At the interface between these domains, the higher-order terms in 
( |lo| ) are important and there is an internal transition layer (similar to a boundary layer). These transition layers can 
be studied using matched asymptotic expansions as has been done for the Cahn-Hilliard equation and other related 
equations similar to (0) jj^] . 

Equation (|l^) is only valid as long as the gradients in the wavenumber are small. During phase slips this is not 
the case. In fact, the phase <\> is n °t even defined at one instant during the phase slip process. The solution shown 
in Figure |l3| includes phase slips and thus cannot be described by (|l0|). However, because each phase slip in the 
solution in Figure |l^ is immediately followed by a second phase slip which restores the phase to its value before the 
first phase slip, on a long time scale the effective total phase A<fi — qdX of the system is conserved. Here the 
effective wavenumber q is defined to be dx4> where </> is the phase of a running time average of A{X, T): 

1 r AT 

A(X, T) = — A(X, T)dT = Re**. (11) 

^ J Jt-at 



The time interval AT is chosen so that flll] ) averages over many double phase slips. An equivalent analysis could also 
be based on the phase of a running time average of B(X, T) instead of A(X, T). The A(X, T) portion of the solution 
to (Q,0) can be written as 

A(X, T) = A{X, T) + Fa{X, T) (12) 

where Fa(X, T) = A(X, T) — A(X, T) contains the rapid fluctuations due to the phase slips and has small magnitude 
except in the neighborhood of a phase slip. A(X,T) contains variations on slower time and space scales and we 



1 Certain phenomena such as the locking of adjacent domain walls into each other are not captured by the phase equation 
([lo|). For small amplitudes of the pattern the locking can be described by suitable Ginzburg-Landau equations |42|^ |. 
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assume that R = R(q) (this is similar to (j6J) ) . Due to the translation symmetry in time and translation and reflection 
symmetries in space of (0,^), the long- wave evolution equation for the effective phase cj>(X,T) is expected to be a 
diffusion equation similar to (jut) 



d T j>(X, T) = D(q)d x 4>(X, T) + h.o.t. 



(13) 



In order to test this explicitly, we measure the effective phase diffusion coefficient D by measuring the response of 
the effective phase of the chaotic state to a localized time-periodic forcing. This approach to measuring the effective 
diffusion coefficient of a chaotic state is similar to that used in experiments on turbulent Taylor vortex flow wherein 
an endcap was moved sinusoidally fig] . Perhaps it should be emphasised that, while we want to measure the effective 
diffusion coefficient to understand the localized chaotic states (as shown in Figure 13), the actual measurements 
described in this section are done on extended chaotic states (as shown in Figure 12). Localized forcing of the 
effective phase in our simulations is accomplished by including a spatially local, time-periodic advection term in (QJ^) 
giving: 



d T A +(v + s)d x A = dd\A + aA + bB + cA{\A\ 2 

d T B + (v- s)d x B = d*d x B + a*B + bA + c*B(\A\ 2 -+ 

v sm(oj T) for X x < X < X 2 
otherwise. 



B\ 2 ) 



gA\B\ 2 

\B\ 2 )+g*B\Af 



(14) 
(15) 



with v 



For small values of vq and ujq, the localized advection term causes formerly stationary solutions to drift within the 
region X\ < X < X 2 , periodically reversing their direction according to the sign of v. For the regions X < Xi, and 
X > X 2 the situation resembles an imposed boundary condition at X = X\ 1 and X = X2 at which the phase of the 
solution varies sinusoidally. 

The solution to (|l3| ) for X > X2, which satisfies the boundary condition (j>(X2,T) = e siu(ujqT), is 



4>{X, T) =Q {X- X 2 ) + ee- a ^ x - x ^ S m(cu T - (3(X - X 2 j) 
2D 



where a = = 

If the effective phase of the solution A(X, T) in the region X > X 2 obeys this relationship we have 



where <f> = ee - Q (*--^) sin(a , oT _ p(x - X 2 )) 
and R = R(d x 4>(X> T )) = &(Qo + d X (j>). 



(16) 
(17) 



(18) 
(19) 
(20) 



If vq and wo are small then the magnitude of the periodic variations in the effective phase 4> will also be small. Taylor 
expanding R(Qo + dx<j>) and e 1 ^ gives 



A{X,T) ~ {R{Q ) + R'{Qo)dx4> 
Taking the Fourier transform of A(X, T) in time, 



JQo(X~X 2 ) 



(l + i4> + ...). 



A(X,Q) = - 



JUT 



A(X,T)dT = - 



AQ.T 



A(X, T)dT + - / e lin F A (X, T)dT, 



(21) 



(22) 



and assuming that the fast variations Fa(X,T) give only negligible contributions to the Fourier components at low 
frequencies yields 



A(X, -Lu a ) 
A(X,lo q ) 



(R(Qo)-0R'(Qo)) + iaR'{Q o ) 
{-R(Qo) - PR'(Q )) - iaR\Q 



,-a(X-X 2 )+i(Q -P)(X-X 2 ) 



-a(X-X2)+i(Qo+P)(X-X 2 ) 



(23) 
(24) 



Thus, the magnitude of the Fourier modes A(X, —u>q) and A(X , luq) (corresponding to the frequency loq in the localized 
advection term) will decay exponentially in space as X — X 2 increases if the phase behaves diffusively. Figure [l6| 
shows that this exponential behavior is realized in our simulations. It gives the absolute value of the Fourier mode 
A(X, wq) versus space for two different values of the frequency u>q. According to ([l7]) the decay rate a is proportional 
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to lUq^ 2 for diffusive behavior. The ratio of the slopes of the curves in Figure are consistent with this diffusive 
scaling. Note that for larger frequencies, and therefore larger gradients of <fi, higher order terms (like —Gd\4>, etc.) 
have to be kept in the diffusion equation which affect the dependence on uj. For G > 0, a is expected to go to zero 
slower than oj 1 / 2 : 

a{uj 2 ) \uj2j 2D 2 

We found analogous behavior in the nonchaotic regime where our simulations were able to reproduce the analytical 
value of the diffusion coefficient. 

The phase slips introduce noise which can be seen in Figure [16[ In order to reduce the effect of this noise and 
extract a reliable decay rate a, the Fourier integral ( p2|) has been extended over many periods (r = 88 x 2it/uiq) of 
the function v(X,T) in the advection term in (14 lj) 



From the decay rate a we can calculate the diffusion coefficient D using (|l7|). We can also calculate D using j3 
which is obtained from the ratio of the Fourier coefficients A(X, u>o) and A(X, —luq) given by (p3|p^): 



& \A(X,-uj q )J 5 I (R(Q )-pR'(Q )) + iaR'(Qo) ' 



A plot of the imaginary part of this quantity versus X will be linear with slope 2/3 if the effective phase is behaving 
diffusively. Figure [l7] shows this plot for the same simulations as presented in Figure [l6| and exhibits the expected 
linear behavior near X\ and Xi- 

The diffusion equation (|l3|) yields a = (3 (cf. (|l7j)) so the slopes indicated on Figure [l7] should be twice those 
indicated on Figure [l6| Although the small discrepancies could be due to higher order terms in (|l^) they are within 
the margin of error in the measurement which is mainly due to the short-scale chaotic dynamics. 



We conclude from the results in Figures |16| and 17 that the dynamics of the effective phase are indeed diffusive 
and show in Figure [l8] the wavenumber dependence of the effective phase diffusion coefficient D(q) along with the 
analytical phase diffusion coefficient for the stationary solution |35|. The solid curve shows the analytical phase 
diffusion coefficient D for the stationary solution. The triangles give the effective phase diffusion coefficient D for the 
chaotic solutions as measured by using the decay rate a, while the circles show D as measured using /?. Despite the 
scatter in the data it is clear that D(q) decreases with decreasing q and presumably goes to zero at a wavenumber for 
which the diffusion coefficient of the nonchaotic state is still negative as shown by the dashed line. Thus the phase 
of the system is diffusively unstable over a finite range of wavenumbers. This, along with the conservation of the 
total phase (here the conservation of the total effective phase), is known to lead to stable domains of differing wave 
numbers, and explains the localization of the chaotic region in Figure |l^. 

For the wavenumbers at which we were able to measure the effective diffusion coefficient, the effective phase of 
the solution remains almost stationary in the absence of any localized forcing. For larger wavenumbers, however, 
the effective phase drifts significantly on the time scale over which we attempted to make our diffusion coefficient 
measurements. Because of this we were unable to measure the effective diffusion coefficient for larger wavenumbers 
than those included in Figure [l^. 



Figure 18 suggests that the chaotic activity should not disappear homogeneously if the average wavenumber is 
decreased towards the stable regime. Instead, the homogeneously chaotic state should become unstable to long-wave 
modulations and eventually split up into chaotic and stationary domains in which the local wavenumber is in the 
respective stable regimes (D > and D > 0). A similar breaking up of the homogeneous state is expected as the 
forcing amplitude b is increased. Numerical simulations confirm this as shown in Figure |l^. 



VI. THE TWO DIMENSIONAL PROBLEM 

In this section we consider a two-dimensional version of (0J|) 

d T A + sd x A = dd\A + d 2 d\A + aA + bB + cA(\A\ 2 + \B\ 2 ) + g A\B\ 2 (27) 
d T B - sd x B = d*d 2 x B + d* 2 d\B + a*B + bA + c*B(\A\ 2 + \B\ 2 ) + q*B\A\ 2 (28) 

These equations describe parametrically excited waves in a two-dimensional anisotropic system. The usual Faraday 
experiment is not described by these equations since it is isotropic. An experimental system that is described by 
(p7],p8|) is electroconvection in nematic liquid crystals [El] , in which the anisotropy is due to a preferred direction of 
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the rod- like molecules. Parametrically excited waves have been observed in both the normal roll regime El , in which 
the waves are perpendicular to the preferred direction of the molecules, and in the oblique roll regime J22], in which 
the waves are at an oblique angle to the preferred direction. The normal rolls are described by ( [27[|2g| ) while the 
oblique rolls are described by four coupled amplitude equations 



The linear stability analysis of the standing wave solution of (^7]28) must consider perturbations transverse to the 
waves as well as along them. This allows additional instabilities as compared to the one-dimensional case. For the 
parameters chosen for the one-dimensional simulations above, along with the dispersion coefficient in the y-direction 
d-2 = d, a stable wavenumber band exists for 0.05 < b < 0.85. 



As in the one-dimensional case (27 28j) can be reduced to a single Ginzburg-Landau equation with real coefficients 
near onset. Therefore the solutions for small b are similar to those of a single real Ginzburg-Landau equation and 
the solution with an initial condition slightly outside the Eckhaus stability limit undergoes a transition which reduces 
the wavenumber. In two dimensions the transition from an Eckhaus unstable wave number to a stable wave number 
involves the creation of defect pairs which are created together and then move apart as shown in Figure ^0|. There, 
part of the spatial domain of the problem at eight sequential times is represented by the eight squares. Within each 
square the solid curves represent locations where the real part of the solution is zero, and the dashed curves represent 
locations where the imaginary part of the solution is zero. Locations where a solid curve crosses a dashed curve are 
the "defects" where the amplitude of the solution is zero. These defects are indicated by circles in Figure |2(]. The 
defects, which form as a pair (at some 500 < t < 600 in Figure p(i| ), move apart leaving one less wave in the domain 
after they have moved far apart (t = 2500). This is analogous to a phase slip in one dimension. 

For larger values of b the solutions to ( p7| , p8| ) also involve the creation of defect pairs. However, the defects often 
do not move far apart but stop after separating some finite distance, turn around and annihilate each other. This can 



be visualized by considering only the left column of Figure 20 , The solution begins with the downward sequence of 
squares from the top as sequential in time (as in the standard Eckhaus transition) . However, on reaching the bottom 
of the column (f = 650) the time sequence reverses and leads back to the top of the column where the wavenumber 
resumes its original value. Such a bound defect pair is analogous to a double phase slip in one dimension. In our 
two-dimensional simulations we have not found the single bound defect pair or simple periodic solutions that are 
suggested by this visualization. But we have seen solutions where many of the defects which are created together as 
pairs annihilate each other. Figure [2l] is a space-time diagram of such a solution for b = 1.2. The spatial aspect of 



Figure 21 is a projection; the y coordinate is not represented. Each circle represents the x-coordinate of the position 
of a defect and the time it was at this location. Small and large circles denote defects of opposite polarity (as they 
do in Figure . Focusing on the "bubble" formed by the pair of defects created at x w 29 and t w 2680 we see that 
the defects in this pair indeed move apart for a short time before they reverse their directions, ultimately annihilating 
each other at x ~ 30 and t rs 2710. There are many such bubbles on the diagram as well as more complex events such 
as bubbles involving two or three defect pairs. While for this value of b most defects are "bound" in such a bubble, 
there are also a few isolated "unbound" defects as well. In these preliminary calculations we have not investigated 
the statistics of bound and unbound defects and the dependence of the statistics on the forcing amplitude b. 



Figures 20 and W] are based on A(X,Y,T). Similar to the phase slips in the one-dimensional problem, in our 



two-dimensional simulations a defect in A is always accompanied by a defect in B at a nearby location. 

VII. CONCLUSIONS AND DISCUSSION 

We have presented numerical solutions of a pair of coupled Ginzburg-Landau equations that are amplitude equations 
for parametrically excited waves ( (p]|2]) , (p7 ,p8f) ) . In particular, we have studied the ensuing dynamics when solutions 
of these equations corresponding to standing waves which are phase locked to the external forcing are perturbed. 
Parameters in the equations were chosen such that there is a forcing amplitude above which all wave numbers are 
linearly unstable (Figure Q). The usual transition from an Eckhaus-unstable state to a stable state (Figures ||a, |2p| ) 
is found near onset but not for larger forcing amplitudes. 

In one dimension double phase slips (Figure [2Jd) occur with these larger forcing amplitudes. These double phase slips 
allow periodic dynamics with an average wavenumber that is Eckhaus unstable even for parameters for which there 
is a range of stable wavenumbers (Figures H g-^Oj). Solutions with an average wave number that is stable can also 
exhibit transient and periodic solutions that include double phase slips. In this regime the system can be excitable 
(Figure |). 

Also found in one dimension are solutions where the dynamics are chaotic and characterized by the irregular 
occurrence of double phase slips in both space and time (Figures [H], |l2| ) . Particularly striking are solutions in which 
the chaos is confined to only a part of the homogeneous system (Figures |l3 - 15 ) . An effective phase diffusion equation 



utilizing the long-term phase conservation of the solution explains the localization of this new form of amplitude chaos 
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(Figures 

A double phase slip in one dimension corresponds to a bound defect pair in two dimensions. A bound defect pair 
consists of two dislocations which are formed together, move apart, then move back together and annihilate each 
other. These bound defect pairs are prominent in our two-dimensional simulations (Figure |2l|). 

So far we have not identified a particular physical system which exhibits the behavior described here. It is natural 
to suspect that the closing of the Eckhaus curve (cf. Figure 0) is related to the occurrence of double phase slips. If 
that should be the case, parametric driving of waves which are Benjamin-Feir unstable in the absence of any forcing 
should be a good candidate for this behavior. Within the framework of the coupled Ginzburg-Landau equations (1,2) 
these waves become unstable near the band center when forced sufficiently strongly |48| . Depending on the sign of 
certain nonlinear coefficients this indicates instability at all wavenumbers and a closing of the Eckhaus curve. Further 
work is needed to investigate the origin of the double phase slips and parameter values for which the dynamics seen in 
our study can be observed. Of particular interest is whether double phase slips might occur in a Faraday experiment. 

Independent of experimental realizability, a very interesting theoretical question arises with regard to the bound 
defect pairs in two dimensions. In thermodynamic equilibrium the unbinding transition of defect pairs has been 
studied in great detail in the context of two-dimensional melting In the vicinity of that phase transition 

quite unusual scaling properties were found. It has been conjectured that such a transition should also be possible in 
nonequilibrium systems [pl"l , [52| |l|l but no very convincing instance has been identified so far. 

The use of an effective phase diffusion equation for the large scale behavior of the chaotic state where the effective 
diffusion coefficient is positive but the usual phase diffusion coefficient is negative (indicating that the non-chaotic 
solution is unstable) is similar to studies using the Kardar-Parisi-Zhang equation [ p3| , which is Burger's equation with 
a Gaussian forcing term added, to model chaotic solutions of the Kuramoto-Sivashinsky equation |54-pq|. Presumably 
a similar noise term £(Y, T) could be added to ( p"3| ) to describe the effect of the shortscale chaotic dynamics on the 
effective phase, 

d T 4> = D d%4> + Didxfa&t + ({X, T) + (29) 

The form of the nonlinear term, which results from expanding the diffusion coefficient D around some wavenumber qo, 
differs from that found in the Kardar-Parisi-Zhang equation since the parametrically driven standing waves have an 
additional reflection symmetry as compared to the traveling waves described by the Kuramoto-Sivashinsky equation. 



Power counting |57| indicates that in the limit of very long scales the nonlinear terms in (J29j) are not relevant (neither 
in one nor two dimensions), in contrast to the nonlinear terms in the Kardar-Parisi-Zhang equation. A particularly 
interesting feature of (|2^) is, however, that even the effective diffusion coefficient Dq can become negative. This is 
not reported for the Kardar-Parisi-Zhang equation as applied to the Kuramoto-Sivashinsky equation. It suggests that 
noise may be quite important in the regime in which D is small, i.e. near the onset of domain formation. 

We would like to thank J. Eggers, G. Grinstein, and Y. Tu for stimulating discussions. This work was supported 
by the United States Department of Energy through grant DE-FG02-92ER14303. 
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FIG. 1. Linear stability diagram for a = —0.05, c = — 1 + 4i, d = 1 + 0.5i, s = 0.2, g — — 1 — 12i. 



Single phase slip Double phase slip 




U=Real((A(XJ)+B(X,T))e iqx ) 



FIG. 2. Space-time diagram illustrating (a) single phase slip observed for small values of the forcing amplitude (e.g. b = 0.1), 
and (b) double phase slip observed for larger values of the forcing amplitude (e.g. b = 0.6). 
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FIG. 3. Example values of the forcing amplitude b and wavenumber Q where the three solution types occur. 




FIG. 4. Phase plane representation of a transient solution involving one double phase slip. The square indicates the initial 
condition which is the perturbed 5-wave state. The loop is traversed once, with the final state of the solution being the (stable) 
fixed point indicated by a diamond, (b = 0.54, Q — 5(2tt)/L — 0.3628). 
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FIG. 5. Each circle represents the location (in space and time) of a double phase slip for a simple periodic solution (6 = 0.54, 
Q = 5(2tt)/L = 0.3645). 
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FIG. 6. Spatial average of the local wavenumber plotted as a function of time for a simple periodic solution (b = 0.54, 
Q = 5(2tt)/L = 0.3645). 
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FIG. 7. Phase plane representation of a simple periodic solution (b = 0.54, Q — 5(2n)/L — 0.3645). 
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FIG. 9. Location of double phase slips for a periodic solution occurring after a period doubling bifurcation (b = 0.54, 
Q = 5(2tt)/L = 0.3653). 
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FIG. 10. Location of double phase slips for a solution where a wave of double phase slips propagates through the system 
(b = 0.54, Q = 5(2tt)/L = 0.3674). 
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FIG. 11. A complex solution where the double phase slips occur irregularly in space and time (6 = 0.54, 
Q = 5(2tt)/L = .3692). 
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FIG. 12. Extended spatio-temporal chaos in a large system (6 = 0.6, L = 800): (a) Q = 50(2tt)/L = 0.3927, (b) 
Q — 51(27r)/L = 0.4006, (c) Initial wavenumber Qo = 52(27r)/L = 0.4084, after a single phase slip at t » 13000 the wavenumber 
is Q = 0.4006. 
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FIG. 13. Localized spatio-temporal chaos in a large domain. Each dot represents a double phase slip (6 
Q = 47(2tt)/L = 0.377). 
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FIG. 14. Temporal evolution of the width of the localized region exhibiting spatio-temporal chaos for two different initial 
conditions (6 = 0.605, Q = 47(2tt)/L = 0.377). 
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FIG. 15. Local wavenumber of the time average of A(X, T) for a localized chaotic solution. The local wavenumber in the 
chaotic region is larger than in the quiescent region. (6 = 0.6, L — 783, Q = 47(2iv)/L = 0.377). 
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FIG. 16. Decay of the magnitude of the coefficient of the time Fourier transform of A corresponding to the frequency ujq in 
the localized advection term v in equations ( |l^Jl5| ), for two different values of uo (b ~ 0.61, L = 800, Q = 51(2n)/L = 0.4006, 
vq = 1.630wo). 
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FIG. 17. Imaginary part of the ratio of the coefficients of the time Fourier transform of A corresponding to the frequencies ±ti; 
in the localized advection term v in equations ( ^4|Jl^ ), for two different values of u>o (b = 0.61, L = 800, Q = 51(2tt)/L = 0.4006, 
vq = 1.630wo). 




FIG. 18. Analytical phase diffusion coefficient D for the stationary solution (solid line) and effective diffusion coefficient D 
for the chaotic solution; Triangles and circles represent numerical results as obtained from a and j3 (cf. Figures |l6| and p^) , the 
dashed line is the presumed extension of these results (b = 0.61). 
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FIG. 19. A set of solutions obtained by incrementally increasing the forcing amplitude b. The extended spatio-temporal 
chaos breaks up into localized regions before the Eckhaus stability limit is reached (L = 800, Q = 51(27r)/L = .4006). (a) 
b = 0.62, (b) b = 0.64, (c) 6 = 0.66. 
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FIG. 20. Standard wavenumber reducing transition in two dimensions. Solid curves represent locations where the real part 
of the solution is zero. Dashed curves represent locations where the imaginary part of the solution is zero. A bound defect pair 
can be visualized by considering the left column only, progressing from the top to the bottom then returning to the top. 
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FIG. 21. Space-time diagram showing the projection of the location of defects on the rr-axis for a simulation of the 
two-dimensional system. The bubbles, for example the one starting at x « 29 and t « 2680, represent bound defect pairs 
(6= 1.2, L = 62.832). 
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